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ABSTRACT 

The herein presented analytical framework fully describes the motion of copla- 
nar systems consisting of a stellar binary and a planet orbiting both stars on 
orbital as well as secular timescales. Perturbations of the Runge-Lenz vectors 
are used to derive short period evolution of the system, while octupole secular 
theory is applied to describe its long term behaviour. A post Newtonian cor¬ 
rection on the stellar orbit is included. The planetary orbit is initially circular 
and the theory developed here assumes that the planetary eccentricity remains 
relatively small (e .2 < 0.2). Our model is tested against results from numerical 
integrations of the full equations of motion and is then applied to investigate the 
dynamical history of some of the circumbinary planetary systems discovered by 
NASA’s Kepler satellite. Our results suggest that the formation history of the 
systems Kepler-34 and Kepler-413 has most likely been different from the one of 
Kepler-16, Kepler-35, Kepler-38 and Kepler-64, since the observed planetary ec¬ 
centricities for those systems are not compatible with the assumption of initially 
circular orbits. 

Subject headings: Celestial mechanics, planets and satellites: dynamical evolution and 
stability, binaries: general 
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INTRODUCTION 


Many problems in dynamical astronomy and celestial mechanics can be reduced to 


describing gravitational 
asteroid (e.g 


interactions in hierarchical triple 


star configurations (e.g 
as black holes (e.g. 


system s. Exam pl es ran ge from 



Borkovits et al. 

Blaes et al. 

2002 

)• 


ary 


systems (e.g. 


Te yssand ier et al. 


2013 1 to multi 


200711 or even systems with compact objects such 


2002). Hierarchical triples consist of a binary system and 


a third body on a wider orbit. The dynamical behavior of such a configuration can be 
pictured as the motion of two binaries: the two inner bodies (inner binary) move on a 
perturbed Keplerian orbit, whereas the third body with the centre of mass of the inner 
binary form the so-called outer binary. Both orbits exhibit coupled dynamical evolution. 


The dynamical behaviour of hierarchical systems, especially the 


behav i our, h a s been a subj e ct of intense study in th e past. 


Kozail (1196211 , iLidovl (1196211 , iHarringtonl (119681 1 and iHeppenheimerl (ll978ll are excellent 


long 


term (secular) 


Brouwer! ( 1959 1. Kaulal ( 1962 1. 


examples of that. Research on the subject has c ontinued throughou t the years with 


Krvmolowski fe Mazeh ( 199911 


( 2010 ), 


Naoz et al 


(l2013al 


Jd), 


Ford et al 


Leung & Lee ( 2013 1. 


2000 ), 


Lee fe Pea. 


Li et al 


e (200311. 


Farago fe Laskar 


(20141) and 


Liu et al. 


(2Q15) 


being a small sample of recent developments in the subject. 

If the system to be modeled is coplanar, then the dynamical evolution of the 
hierarchical triple is dominated by the time dependent changes in the eccentricity vectors 
(Laplace-Runge-Lenz vectors) of the inn er an d outer orbit, as the semi-major axes do 


not show any secular evolution (e.g. 


Marchal 


199011 . Hence, a solution of the differential 


equations governing the behaviour of the eccentricity vectors allows us to accurately 
describe the motion of the whole system via closed analytic formulae. 

In a series of papers, we f ocused o n the evolution of the inner eccentricity of a 


hierarchical triple system. In 


Georgakarakosl ( 20021) . we developed a method for measuring 
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the eccentricity injected to the inner orbit due to the perturb ations of the third bo c 


effect has been known for years (e.g. 


Mazeh fe Shaham 


1979), but 


Georgakarakosl (120021 ) 


y. This 


quantified it analytically on both short term and secular timescales for the first time. 
As those calculations were performed for coplanar and initially circular or bits, t he y were 


l ater on extended to cases where the perturber was on an eccentr ic orbit (IGeorgakarakos 


2QQ3|) as well as for three dimensional cases (IGeorgakarakos 20041) . The initial motivation 
for th at work was close ly related to modeling tidal friction in stellar triple systems (e.g. 


see 


Kiseleva et ah 


1998!) . Therefore, the systems investigated had comparable masses. 


Nonetheless, the derived formul ae wer e _also tested numerically for planetary mass ratios 

.iinates in the case of 


and remained valid (IGeorgakarakos 


2006). Finally , an im proved e s 


coplanar and initially circular orbits were given in 


Georgakarakosl (2009). The analytical 


formulae 

versatile. 


'or th e eccent ricity evolution in hierarchical triples have proven to be quite 


Chavez et al.l (20121) used them to investigate the long term modulation in the 


lightcurve of the cataclysmic variable FS Aurigae, which was suspected to happen due to 
the gravitational perturbation of a third body on a much wider orbit. More recently, the 
limits of various types of habitable zones for a planet in a stellar binary with di fferent 


spect ral type components were investigated analytically as well as numerically (lEggl et ah 


2012 


2013). 


Boostec 

by 

le discovery of several circumbinary planets fe.g. 

Dovlc et ah 

2011 


Welsh et ah 

2012: 

Orosz et al. 

2012: 

Schwamb et ah 

2013 

Rostov et ah 

2014), the interest 


in analytical models that describe cicumbinary planetary motion has recently grown. In 
order to create more accurate analytical tools for detection algorithms and planetary 
formation theories we have to improve first and foremost our description of the eccentricity 
evolution of planets on circumbi n ary o rbits (also known as P-type orbits). There has been 


some progress m 


Georgakarakosl ( 20091 ) regarding the latter matter, but only hierarchical 


triple systems on initially circular orbits were treated therein, which may be inadequate in 







































































5 


certain c ases. The Kep l er-3 4 system for example has a stellar binary with an eccentricity of 


e ~ 0.52 flWelsh et al 


201J). 


In this paper, we derive analytical expressions for the eccentricity of a relatively light 
body orbiting a heavier binary, e.g. a planet moving around an eccentric double star. 
This will allow us to model the dynamical evolution of a circumbinary planet analytically. 
Since most of the circumbinary planets that have been found via transit searches share 
more or less the same orbital plane, we will focus on coplanar configurations. Furthermore, 
we include a post-Newtonian correction into our model, since, in certain cases, General 
Relativity (GR) may affect the orbit of the inner binary, which in turn influences the 
evolution of the planetary eccentricity. Our model assumes that the planetary eccentricity 
is initially circular and it remains valid for low values of the planet’s eccentricity (e 2 < 0.2). 

The remainder of this article is structured as follows: in section [2 we present the 
analytical derivation for the eccentricity of the planetary orbit, as well as the expressions 
for its maximum and averaged value. In section [3], the analytical estimates are compared 
to results obtained from the numerical integration of the full equations of motion. Next, 
we discuss the effect of a post-Newtonian correction in our problem (section [4| and we 
present the formulae for the analytic propagation of the hierarchical triple system (section 
[5]). In section [6l we apply our analytical estimates on six circumbinary planetary systems 
discovered during the Kepler mission. Finally, we discuss and summarize our results in 
section 0 A brief description of the variables used in this work can be found in the notation 
section after the appendix. 
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2. THEORY 


We are interested in finding a complete description of the eccentricity vectors of the 
inner and outer orbit, since those will determine the dynamical evolution of a coplanar 
hierarchical triple system. As it was stated earlier, we will assume a star-star planet 
configuration, where all bodies are treated as point-masses. Since the outer body is of 
planetary size, we can neglect its influence on the amplitude of orbital eccentricity of 
the inner pair (this will be justified in section |3]). Only changes in the direction of the 
eccentricity vector of the inner orbit will be taken into consideration. 


The general way to derive accurate eccentricit y est imates that contain long and short 


periodic contributions has been laid out in 


Georgakarakosl (120031) . First, the long term 


evolution of the eccentricity vector (Laplace-Runge-Lenz vector) is studied by deriving and 
solving equations of motion that have been averaged over all fast angles. The averaging 
process via a canonical transformation eliminates all short periodic terms in the eccentricity 
evolution, so that they do not appear in the final solution. If we want to describe the 
system’s behavior on orbital timescales as well, we can reintroduce short periodic terms by 
deriving the non-averaged equations for the eccentricity vector. We then eliminate all the 
components that are already contained in the long term (secular) solution and solve the 
reduced differential equations. The final solution is a superposition of the short and long 
periodic solutions of the eccentricity vector. 


That kind of approach has been chosen in this work. The secular solutions of the 
outer eccentricity vector have been constructed from the equations of motion derived from 
the corresponding canonically averaged Hamiltonians. The short periodic terms, however, 
had to be themselves split into medium (with periods of the order of the planetary orbital 
period) and short term (with periods of the order of the stellar orbital period) components 
in order to allow for an analytic solution. Rejoining the different contributions of the 
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long, medium and short term solutions leads to the final result as detailed in the following 
sections. 


2.1. Calculation of the long term evolution 


In order to de rive the long-term modulation of the system, we shall use the Hamiltonian 


formulation ( Marchal 


19901 ). For this purpose, we make use of the Delaunay variables, 


which is a set of canonical variables. For a coplanar three body system they are defined as 
follows: 


( 1 ) 


mnia\ , 

h, 

Ain 2 a 2 , 

h, 


9i = tti, 

L 2 \/l - e\, 

g 2 = zu 2 , 


where a*, e,, tu* and are the semi major axes, the eccentricities, the longitudes of 
pericenter and mean anomalies of the inner (i = 1) and the outer (i = 2) orbit respectively. 
Furthermore, 


m = 


m 0 mi 


and A4 = 


m 2 (m 0 + mi) 


m 0 + m 1 M 

are the so-called reduced masses, with mQ,mi and m 2 , being the masses of the stellar 

binary and planet respectively. The total mass of the system is M = Y2i=o m i- Using those 
variables, the Hamiltonian for a hierarchical triple system reads 


H = H 0 + H x + H p , 

where 

_ Q 2 mlml 

0 2(m 0 + mi)L^ 

is the Keplerian energy of the inner orbit, 

^ _ G 2 (m 0 + mi) 3 ml 

1 ~ 2 ML 2 


( 2 ) 

( 3 ) 

( 4 ) 











is the Keplerian energy of the outer orbit, and 


rj n ,TO 0 + mi m 0 m 1 

H p = gm 2 { ---), 

R r 02 r 12 


(5) 


is the perturbing Hamiltonian, with r 02 and r 12 being the distances between m 0 and m 2 and 
mi and m 2 respectively. Furthermore, R is the vector which connects the centre of mass of 
the inner binary with the third body, while r will be denoting the relative position vector 
of the inner orbit (see Figured]). Finally, Q is the gravitational constant. 

Since we deal with a hierarchical system and r/R is small, H p can be expressed in 
terms of Legendre polynomials, giving: 

Qm 0 mim 2 


n„ - - 


R 


£v,(-)^(cos9), 


3 =2 


( 6 ) 


where Mj = —, Rj are the Legendre polynomials and 9, as seen in Figured! is 

the angle between the r and R vectors. 

The long term behaviour of the system can be obtained when the sh ort pe riod effects 


in the Hamiltonian are removed by using the Von Zeipel method (e.g. 


Marchal 


19901 ). The 


method e mployes a canonical trans 


elsewhere flBrouwer 


1959 


Marchal 


'ormation and details about t 


1990 


Krvmolowski & Maze 


d 


re derivation can be found 


1999! ). The advantage of 


the Von Zeipel transformation compared to other averaging methods, such as for example 
the ’scissors method’, is that it produces a Hamiltonian which is second order in terms 
of the masses (although that term will not be used here as its effect is negligible due to 
the fact that the third body has a planetary mass) and which is complete in terms of the 
eccentricities. Hence, the Hamiltonian for coplanar orbits averaged over the inner and outer 
mean anomalies is 


< H >—< H 0 > + < Hi > + < H p2 > + < H p 3 > +0 (—), 


(7) 
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with 


< H 0 > = 
<H 1 > = 


Q 2 m^m\ 


2 (m 0 + mi)L^ 
Q 2 {rriQ + mi) 3 ml 
2 ML| 


rr ^ _ l G 2 (m 0 + m 1 ) 7 ml L\ t o G\^ 

P2 § mgmfAf 3 L|G| L? { 

15 G 2 {rriQ + mi) 9 m2(mo — mi) L b v 


<H P 3 > = 


^ 21^*21 


/~i 2 
Lt 11 

L 2 u 


r <2 

( - JC 2 l 

Lit 


X 


( 8 ) 


64 m^m^M 4 

G 2 

x (7 - 3—7p) cos (^i/ - # 2 z) 

L u 

where the index l stands for the long term (secular) evolution and H p2 and and //p 3 come 
from the P 2 and V 3 Legendre polynomials respectively. 

Using Delaunay variables, Hamilton’s equations assume the following form for a 
coplanar system: 


dLji d < H > 

dt dlji 

dGji _ d<H> 

dt dgji ’ 


dlji d < H > 

dt dLji 

dgji d < H > 

dt dGji 


( 9 ) 


Since our Hamiltonian is independen t of lu- all Lu ar e constant. This is a well known 
result for the averaged probelm (e.g. 


Ha rr ington 


19681 ). Then, keeping in mind that the 


eccentricity of the stellar orbit remains unchanged due to the fact that the outer body has 


















a planetary mass, the relevant equations of motion are: 


dG^i 

dt 


dg 2 i 

dt 


dgu 

dt 


d<H p3 > 

dg2i 


15 G 2 (rn 0 + — mi) L b h 


u 


64 


momfM 4 





x - 


' 1 - jf (7 - 3jt) sin (dii ~ gn) 
^21 ^ '' 


d < II, 


P 2 


Hp 3 > 


3 G 2 (m 0 + mi) 7 m 7 2 


d G 2 i 8 rriQm\M' i 

15 G 2 (m 0 + mi) 9 ml(mo — mi) 

64 


T 4 n 2 

r3 /^<4 0 t 2 > 

^ 21^21 ^11 


m^m\M 4 




/^2 rr 2 

^(7-3^) 


X 


J ll 


J ll 


X 


'5-i/l 


G|, 


G 2l 


LhG&J 1 


"Gi: 


cos ( gu - g 2 i) 


d < II, 


p2 


H P 3 > 


3 G 2 (m 0 + mi) 7 ml L\ t Gn 


dGu 4 m^mlM 3 

15 G 2 (m 0 + mi) 9 ml(m 0 — mi) 


L~2l G 2l 


+ 


m2 
^21 


G 1 


m^m^M 4 


64 

G‘2 

■9-^) cos (^H - g 2 {) 

L \i 


T 3 /^5 
lj 2l Ksr 2l 


L 2iv'W zr G I u 


(13 


( 10 ) 


( 11 ) 


( 12 ) 


If we replace the Delaunay variables from equations (JTJ) with the corresponding orbital 
elements and use e 2 u = e 2 icosg 2 i and e 22 i = &21 sin g 2 ;, the above system becomes: 
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de 


21 1 


dt 


de 


22 1 


dt 


dgu 

dt 


y/QMmomiah 


2 12 
21 ) 


X — 

° (m 0 + mi) 2 a 2l (l — e 
15 \/GMm 0 mi(m 0 — mi)a 3 z 


-(2 + 3e\i)e22i + 


64 


(m 0 + mi) 3 a|;(l — e 
+5e 2 i Z e 22Z cos^ij 
3 V QMmomiah 

2 12 


2 13 
2/1 


e iZ (4 + 34) [(1 + 4e2 2Z - e 2 21l ) sin g u + 


(mo + mi) 2 a 2 Z (l — e; 

15 \/QMmomi(m 0 — mi)a 3 


-(2 + 3e^)e 2 i Z + 


64 


(m 0 + m 1 ) 3 a| i (l - e|) 3 


—5e 2 i Z e 22Z sin^i;] 

3 

3 VQm 2 al l y/T 


- 1 1 


_15 VUm 2 (m 0 - m^a^ea^/l - e 2 , 

4 (m 0 + m!) 2 a|(l - e |)2 64 (m 0 + m^a^euy/1 - e| 

+94) cos((?i Z - g 2 i). 


(4 + 


(13) 


■eiK 4 + 34) [(-1 + e\ 2l - 4e4) cos g u - 


(14) 


(15) 


Since the planetary orbit is initially circular and its eccentricity is expected to remain 
small, we neglect terms of order 0 (e 2l ) in the above equations and we also keep the 
dominant term in the equation for the stellar pericentre. As a result, the equations of 
motion assume the following form: 


where 


de 


21 1 


dt 

de 22 i 

dt 

dgu 


dt 


= -K\e 22 i + K 2 sin g lh 
= K\e 2 u — K 2 cos gu, 

= k 3 , 


3 VGMmom^h ^ 2 

Ai = - r (2 + 3e u ) 

° (m 0 + mi) 2 a 2l 
15 \/QMmomi(m 0 — mi)a' 


K, = 


64 


(m 0 + mi) 3 ^ 


e iZ (4 + 3e 


K, = - 


3 \fQrn 2 a 2 u ^J\ — e 


4 (mo + mi)2a 3 z 


(16) 

(17) 

(18) 


(19) 
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are constants. 


The above system of differential equations can be solved analytically and the solution 
is 

K 2 

e 2 u(t) = C i cos Kit + C 2 sin K 1 t+ — -— cos K 3 t 

Ki — A 3 

k 2 

e 22 i (t) = C i sin K\t - C 2 cos K\t + ——sin K 3 t, 

A i — A 3 

for the planetary orbit, where C\ and C 2 are constants of integration. For the stellar orbit 
we get 

9 u = K :i t, (21) 

where we have assumed without loss of generality that the initial longitude of pericenter of 
the inner orbit is zero. Consequently, the components of the eccentric vector of the stellar 
orbit become 



eiu (t) = ei cos g u = e x cos (K 3 t) 
ei 2 1 (t) = ei sin g u = ei sin (K 3 t ), 


( 22 ) 


So f ar ou r resul ts are in agr eement with other authors, e.g. 


(20041). Demidova fe Shevchenko (2014). 


Moriwaki . & Nakagawa 


2.2. Calculation of the short and mid term evolution 


Having the long term solution for the eccentricity evolution of our system at our 
disposal, we now proceed to derive equations for the mid and short periodic terms. Once 
again, we use the Jacobi decomposition of the three body problem to describe the motion 
of the system (see Figured]). In that context, the equation of motion of the outer body is: 


it 


-QM 



R + fi\r 
R + /iir | 3 


+ /A 


R-^or \ 
R — no n| 3 / 


QM 


d 




ho 


+ 


hi 


dR\\R + Hir\ \R —nor\ 


(23) 
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with 


rrii 

Vi =-—-, 

mo + m i 


i = 0,1. 


Once more, since we work with a hierarchical triple system and the third body is 
assumed to be at considerable distance from the inner binary, implying that r/R is small, 
the inverse distances in equation (|23j) can be expressed in terms of Legendre polynomials as 


1 

n —0 

(24) 

\R + ni r \ 

and 



1 

n =0 

(25) 

\R — hqv\ 


Expanding to second order and making the substitution 


cos 6 = 


r- R 

rR ’ 


the equation of motion of the outer binary becomes 

_ r R r (r- R)r 

R = GM — — + HoVi [3- 


R 5 


15 f r-R) 2 R 3 r 2 R 


R 7 


2 R 5 


(26) 


Using the Runge-Lenz vector of the outer orbit we can obtain an expression for the 
eccentricity. Hence: 


e 2 = 

e 2 = 

where h = R x R. 


R 

R 


QM 


[R x h) 


R + R rt 

~r + 1? r 


QM 


[2 (R R)R - (R R)R - {R R)R], 


(27) 


Substituting equation (]23|) into equation (127|) and assuming that the eccentricity 
remains small (and hence neglecting the term R ■ R), we obtain: 

momi 


^2 sm 


(m 0 + mi ) 2 


r (r-R)(r-R) n 9 (r ■ R) 2 ■ 3 r 2 • n 

[6-- ^ 


ir 


2 R 3 


(28) 
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where the indices s and m refer to short and medium period terms as explained earlier. 

Now, the Jacobi vectors can be represented as r = ci\ (cos E\ — ei, yH — e\ sin £j) and 
R = a 2 (cos / 2 , sin / 2 ), where cq, e\, E\, are the semi-major axis, the eccentricity and the 
eccentric anomaly of the inner orbit, while a 2 and l 2 are the semi major-axis and the mean 
anomaly of the outer orbit respectively. Substituting for r and R in equation (jZ8]) . we 
acquire for the rate of change of the componets of the outer eccentric vector: 


621 sm{t) — 


mouii 


n 2 


4 , ,2 „4 


6 [(cos Ei — ei) cos l 2 + Jl — e\ sin E\ sin l 2 ] x 


(m 0 + mi) 3 M 3 J 3 - 

/- 9 

x [ — (cos Ei — ei) sin l 2 + y 1 — e 2 sin Ei cos l 2 \ cos l 2 -[(cos E x — 


ei) cos l 2 + y 1 — e\ sin E x sin Z 2 ] sin l 2 + -(1 — ei cosi?i) 2 sin l 2 


&22sm (t) 


UIqUIi 


n 2 


4 , r 2 ,,4 


6 [(cos Ei — ei) cos Z 2 + y 1 — ef sin E x sin l 2 ] x 


(m 0 + mi) aM3 X 3 

/- 9 

x [ — (cos Ei — ei) sin l 2 + y 1 — e\ sin E x cos l 2 \ sin l 2 + - [(cos — 
— ei) cos l 2 + y 1 — e\ sin Ei sin l 2 ] cos l 2 — -(1 — ei cos-Ei) 2 cos l 2 


(29) 


where X is the period ratio of the two orbits. 


By using equations (12§|) . we are now able to obtain the short period terms for the 
evolution of the planetary eccentric vector. This will be done in two steps. First, we average 
the above equations over the fast motion and then we integrate over time: 


6 2 1 m(f) — 

&22m(f) = 


1 m^mi 1 

16 ( 777,0 + ■ 

1 m^mi 1 

16 (m 0 + xt 


12 cos l 2 + e 2 (33 cos l 2 + 35 cos 3Z 2 ) 
■ 12sinZ 2 + e 2 (3sinZ 2 + 35sin3Z 2 ) 


+ Ce 2 i m 

F C < e22m . 


(30) 


Here, C e21m and C e22m are constants of integration. In a second step we add higher frequency 
terms by integrating equations (1291) with respect to time. During that process, any terms 
that appear and have already been taken into consideration (given by equations (j30li ) are 
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removed. Integrating by parts yields a power series solution in terms of 1/X. By keeping 
the largest term in that power series we finally end up with the following short term 
equations for the eccentric vector of the planetary orbit: 


e 2 i sit) = 


e 22 sif) — 


rriQmi 


(mo + mi)sM3 Xs 
m 0 mi 1 


■Pi(t) + C e 


(mo + mi) 3Ms X 3 


P2(t) + C, 


e22s 5 


(31) 


with P\ (t) and P 2 (t) given in the Appendix. Combining equations (130|) and (13T)) we hnd 
formulae for all non-secular contributions: 


C 2 1 sm{t) 
f- 22 sm(^) 


mom! 


■ cos l 2 + e\ 


2 /33 , 35 OI . Pi(t) 

2 ' — cos l 2 + — cos 3/ 2 ) H- 


(m 0 + mi)sM3 X 3 L4 1 16 

1 r3 


momx 


(m 0 + mi) 3 M 3 X 3 


7 sin l 2 + e x (— sin l 2 
4 v 16 


16 

35 , , 

— sm 3 / 2 ) + 
16 


X J 

AM 

x 


+ Ce 


(32) 


+ a 


e 22 s 


We would like to point out here that a 1; a 2 and ei have been treated as constants in the 
above calculations. 


2.3. Complete solution of the eccentricities 

We can now combine equations (|20l) and (T32|) in order to obtain the total solution for 
the outer eccentricity. As in the case of the inner eccentricity, this is done by substituting 
the constants of integration in the short period solution with the expressions for the secular 
evolution: 

e 2 i if) — e 2 i s{t) + e2im(t) + e 2 i i(t) — C e21sm 

(33) 

e 22 (t) = e 2 2 s(t) + e 2 2 m{t ) + e 2 2 i(t) — C e22arn . 

The constants in equation (l33l) are determined by the fact that we assume the outer 
eccentricity to be initially zero, i.e. 


e 2 i(to) — e 2 2(to) — 0 


( 34 ) 
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which eventually leads to 
m 0 mi 


e 2 i (t) 


(mo + mi)3Ms X3 L 4 
C\ cos K\t + Ci sin Kit + 


3 , 2/33 , 35 OI v Pi{t) 

cos l 2 + e 1 [— cos 1 2 + — cos 3f 2 ) + 


Ah 


Ah - K, 


cos K 3 t 


e 2 2 (t) 


m 0 mi 


1 


16 


35 


X J 


+ 


u u <ju A 2 (t) 

4-y—4 7 sin / 2 + e? ( 77: sin / 2 + 777 sin 3/ 2 ) H-77- 

(m 0 + mi)lMl Xl U u 16 16 ’ X \ 


+ 


C'i sin Kit — C -2 cos Kit 


K, 


Ki - Ah 


sin AT 3 t, 


where 


Ci = 


momi 


4 , ,2 ,„4 


(m 0 + mJsMs J3 
Ai(to 


■3 , 2 ,33 , 35 . 

7 cos Z 2 o + ej (— cos t 2 o + — cos 3Z 2 o) + 
4 16 16 


Ah 


Ch = 


X 


m 0 mi 


1C - K 3 


(mo + mi)sM3 X 3 L4 

. A 2 (t 0 


3 3 35 

sin ho + (— sin / 20 + — sin 3/ 2 q) + 


16 


X 


(35) 


(36) 


(37) 


The quantities K, have been defined in equations (1191) and Z 20 is the initial planetary mean 
anomaly. The evolution of the eccentricity of the inner orbit is given by equations 


2.4. Maximum and average squared outer eccentricity 

We can also obtain an estimate for the maximum value the outer eccentricity is bound 
to reach during its evolution assuming it was initially zero. This is achieved by doing some 
algebraic manipulations and by maximizing trigonometric functions in equations ( 120 |) and 
(132|) . Eventually we obtain: 

max max , max 

e 2 e 2 sm ' e 2 1 

m 0 mi 1 

7TT~2 _ _ 4 

+ mi) 3M3 X3 

2 Ah 

+ A, - Ah' 



17 , 1 / 21 9 3 n\ 

+ 7 er + v( 3 + 19ei + Y ei "2 e ‘) 


( 38 ) 
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Averaging equations 
short period terms): 


over time we obtain the following expression (excluding the 


( e 2 )t ~ (<22lW+ e 22 (t))t ~ 


2 2 
m^m\ 


75 


— c-l + e\ I 
16 1 1 V 32 


225 


+—el + cos 2/20 T 


(mo + mi) 3 Ms AT3 
525 


9 27 2 887 4 

8 + v' 1 + Hi ' =i + 


128 


cos 4/ 


20 


+ 


m 0 m i 


1 


Ah 


,/33 


(mo + mi) 3M3 A"3 Ai — A 3 

Ah 


2 cos ho + ef y cos Z 20 + 


35 0 / 

COS 3/20 


+ 2 


Ad - Ah 


(39) 


If equations (|35li are averaged over both time and initial angles, we find the most 
compact expression for the averaged square outer eccentricity: 


( 4 ) = <e!i(f) + 4.(*)> = 

975 1 


2 2 
m^mi 


9 27 2 887 4 

(mo + mi)§M§ At |_8 8 64 1 


+2 


64 X 1 

Ah 


1 /225 6619 


X 2 ( 64 


+ 


64 


' e i - 


26309 

"512" 


- e i - 


393 

6 T f 


+ 


Ah - Ah 


(40) 


NUMERICAL TESTING 


In order to test the validity of our analytical estimates, we integrated the full equations 


of motion numerically. For t 
transformation developed by 


l at purpose, w e used the symplectic integrator with time 


Mikkola ( 19971 ). specially designed to integrate hierarchical 


triple systems. The code uses standard Jacobi coordinates, i.e. it calculates the relative 
position and velocity vectors of the inner and outer orbit at every time step. Those were 
used to generate orbital elements of the stellar and planetary orbits. 


We investigated systems with nine different mass combinations. Stellar masses were 
chosen such as mi/(m 0 + m^ = 0.5, 0.3, 0.1. The planetary masses we used were such that 
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^ 2/(^0 + mi) = 10 2 ,10 3 ,10 6 , i.e. we chose Earth-like, Jovian and brown dwarf like 
companions. The eccentricity ranged from 0 to 0.9, while for the semim ajor axis t he range 


was b etween one and 12 times the critical semi-major axis as given in 


Holman & Wiegert 


(119991) . who studied the stability of planetary orbits in binary star systems by performing 


numerical simulations in the framework of the planar elliptic restricted three body problem. 
If a particle survived the whole integration time at all initial longitudes, then the system 
was classified as stable. For circumbinary orbits, and by using a least squares fit to their 
data, they obtained an expression for the critical planetary semimajor axis a c , i.e. the 
smallest planetary semi-major axis for which the particles were stable at all initial positions, 
given by (using our notation for the mass parameters and orbital elements): 

a c = [(1.60 ± 0.04) + (5.10 ±0.05)ei +(-2.22 ±0.U)e 2 +(4.12 ±0.09)^ + 

+ (—4.27 ± 0.17)ei/ii + (-5.09 ± 0.11 )/i 2 + (4.61 ± 0.36)e 2 /i 2 ]ai. (41) 


Finally, the initial stellar eccentric anomaly and the initial planetary mean anomaly 
were varied from 0° to 360° with a step of 45°. The results from the numerical simulations 
were compared with our analytical estimates for (e 2 ) and e™ ax , i.e. equations (1381) and (1401) . 
on both short and secular timescales. For the short timescale simulations the integration 
time was one planetary orbital period, while for the long term simulations the integration 
time was set to 2n/(Ki — ib 3 ), which is the analytical estimate for the secular period of the 
planetary eccentricity. 

Generally, our analytical estimates were in good agreement with the results obtained 
from the numerical simulations. As expected, the greatest errors were for systems where 
the planet was started at the critical semi-major axis. As we move away from the critical 
semi-major axis the error reduces rapidly. Figures [2]{5] show the most relevant numerical 
findings. Finally, Figure [6] is an example that demonstrates the validity of our assumption 
that the semi-major axes and stellar eccentricity remain constant. 
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POST-NEWTONIAN CORRECTION 


For certain stellar binaries the gravitational field between the two stars can be strong 
enough to make it necessary to include General Relativity in order to describe the motion 
of the system correctly. In such cases, our theory can easily be modified to account for the 


extra perturbation. One simply adds the 
the stellar orbit to the equation (J2J) (e.g. 


following first ore 


Naozetah 


2013b ): 


er post-Newtonian correction for 


H\pn — 


(rrip + m\)p A G{3ml + 7momi + 3m\)p 2 G(p-r) 2 G 2 (m 0 + mi) 2 m 


, (42) 


8c 2 mQml 2 c 2 mom 1 r 2 c 2 r 3 2 c 2 r 2 

where p and p are the linear momentum and its magnitude respectively and c is the speed 
of light in vacuum. Averaging the above Hamiltonian over the fast motion, we obtain 


flNaoz et a! 


2013bh : 


< H\pn >— 


+ 29momi + 15m 2 ) 


3 t/ 4 mQ? 7 R 


8(m 0 + mi) 3 c 2 L 4 z (m 0 + mi)c 2 L 3 ; Gi/ 

which eventually results in having the following extra term: 

3 3 

dgiPN 3^2(m 0 + mi)2 


dt 


c 2 al{l 


'ii) 


The above term can simply be added to equation (lT5j) . Then we have 


and consequently 


where 


dgi dgu dgiPN 

- «- 1 - 

dt dt dt 


gi(t) ~ Kzpn t 


Kzpn — - 


3 G 2 m 2 u^ (1 — e 


2 


+ 


3G 2 {m 0 + mi) 2 


(43) 


(44) 


(45) 


(46) 


4 (m 0 + mi)2a 3 z c 2 a^(l - e 2 u ) 

An example of the effect of the post-newtonian (PN) correction is given in Figure [7] 

Here, we have integrated numerically the Einstein-Infeld-Hoffmann equations using a 


( 47 ) 
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Gauss-Radau scheme (jEggl fe Dvorak! 120101 ) for a system with mo = 5M 0 , mi = 4M©, 
m 2 = lMj, ai = 0.2 au, ei = 0.5, tu = 0°, E w = 0° and Z 2 o = 90°. The initial period ratio 
for the above system was X=200, yielding a planetary semi-major axis a 2 = 6.84014558 au. 
The graph demonstrates that the inclusion of GR in modelling the motion of the stellar 
binary can have a significant effect on the evolution of the planetary eccentricity. 

The post-Newtonian correction can be applied to our analytic estimates by simply 
exchanging K :i with K 3 P n, for instance, in equations (135]) and 


5. ANALYTIC DESCRIPTION OF THE SYSTEM’S EVOLUTION 

In this section we demonstrate how equations ( 1221 ) . ( 1551 ) and ( 1451 ) can be used to 
describe the time evolution of the entire system, i.e. the binary and the planet. For the 
inner orbit we have 


w x « gi = K 3PN t, (48) 

h = flit + hoi (4:9) 

rcos/i = ai(cosE 1 — ei), (50) 

r sin /) = ai(l — ei) 1//2 sinEi, (51) 


where fi denotes the corresponding true anomaly. Here we have assumed that di — e± — 0. 
For the mean motion follows that ni = Q 1 ^ 2 ^mo +mi) 1 ^ 2 a 1 3 ^. Defining W- L as the simple 
2D rotation matrix 


W, = 


COS W{ — Sill Wi 

sin Wi cos Wi 


i= 1,2 


we find the solution for the relative vector of the binary stars, r (t) is 


(52) 


r— Wi(r cos fi,r sin fi) T . 


(53) 
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In order to find the position vector of the planet, one follows the same procedure, except 
now the length of the eccentricity vector is no longer constant in time and we have 


w 2 

= arctan(e 22 /e 21 ), 

(54) 

k 

= n 2 t + w 2 + l 2 0 , 

(55) 

e 2 (t) 

= |e 2 | (t), 

(56) 

R cos f 2 

= a 2 (cos E 2 - e 2 ), 

(57) 

R sin f 2 

= a 2 (l - el ) 1/2 smE 2 , 

(58) 


where n 2 = {QM) l ^ 2 a 2 3 ^ 2 . Here, the corresponding values for e 2 have to be taken from 
equations (j551) . The relative vector form the barycenter of the binary star to the planet is 
then given by 

R = w 2 (R COS f 2 ,R sin / 2 ) t (59) 


Again, we have assumed that a 2 = 0. We would like to point out here that the eccentric 
anomaly E % for any of the two orbits can be computed by solving Kepler’s equation 
U = Ei — e t sin E % . Alternatively, a series e xpansio n that rel ates the mean and the true 


anomaly could be used, for instance (e.g see 


Mur ra y & Derm ot t 


19991): 


fi = h + 2 a sin k + ^e 2 sin 2 U + 0(e 3 ) (60) 

However, one should keep in mind that the above series solution is divergent for 
e* > 0.6627434. 


A transformation to barycentric coordinates can be easily achieved through 


r 0 b 

rib 

r 2 b 


m 1 m 2 ,, 

- r - -R 

mo + mi M 

m 0 m 2 

-r—— R 

m 0 + ffi] M 

m 1 + m 2 

-TT- R 


where is the vector of the m* body with respect to the system’s barycentre. 


(61) 

(62) 

( 63 ) 
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6. APPLICATIONS TO REAL SYSTEMS 


We now apply our analytical theory to real circumbinary planetary systems. For that 
purpose, the Kepler-16, Kepler-34, Kepler-35, Kepler-38, Kepler-64 and Kepler-413 systems 
were selected, as they are currently believed to harbour only one planet in a circumbinary 
orbit. The systems are assumed to be coplanar, vj\ = 0°, E w = 0° and Z 2 o = 90°, while 


the rest of t 

le svs 

;em parame 

ers were taken from the corresponding discoverv papers 

(Dovle et al. 

2011; 

Welsh et al. 

2012; 

Orosz et al. 

2012; 

Schwamb et al. 

2013; 

Kostov et al. 


20141 ). The systems were integrated over one analytical secular period and no other effects 
than Newtonian gravity were considered, as they were n ot expected to make a significant 


contribution to the systems under investigation (e.g. 


Chavez et ah 


2015 '). Table 1 gives the 


mass parameters and orbital elements of each system. 

Figures [8] and [9] show the results for the six Kepler systems. Generally, the numerical 
results are in good agreement with the analytical estimates. Furthermore, one can see that 
for most planets the current state of eccentricities, indicated by a black horizontal line, 
is compatible with formation scenarios that predict initial orbits with low eccentricities 
after the gaseous phase. As in-situ planetesimal accretion as well as gravitational collapse 


have practically been ruled out for m ost o f the circum 


Kepler mission (e.g. 


Pelupessy & Portegies Zwa.rt 


j inarv pl ane 


2013 


Is disc overed during the 


Lines et al. 


2 0141 ). a fast disc driven 


migration with little time spent near resonances seems to be the most likely for mat ion 


scena rio for Kepler-16, Kepler-35, Kepler-38 and Kepler-64 (e.g. 
2014j). 


Kiev & Haghighip ou r 


Exceptions are Kepler-34 and Kepler-413, both with a higher planetary eccentricity 
of e 2 = 0.182 and e 2 = 0.1181, respectively. Looking at the relevant plots, it is clear that 
starting the planet on a circular orbit can not produce a planetary orbit with eccentricities 
higher than 0.03 for Kepler-34b and 0.04 for Kepler-413b. Moreover, the main eccentricity 
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contribution for both systems comes from short period activity. This is to be expected, as 
the stellar masses of Kepler-34 have only around 2.5% difference, and the stellar eccentricity 
of the Kepler-413 is just 0.0365. As a result, the forced secular eccentricity, which is 
proportional to the difference between the masses of the stellar components and to the 
stellar eccentricity is very small. Therefore, either those two planets were formed on a 
non-circular orbit or if they were initially circular, some dynamical event may have taken 
place and pumped up their eccentricity. For instance, an as of yet undetected companion 
as well as an encounter with another star may have injected eccentricity into the planet’s 
orbit. Such an interaction would also explain the slight misalignment of the orbital planes in 
Kepler-413. Another possible explanation for the elevated Kepler-34b is resonance trapping. 
If the planet’s migration has not been fast enough, the planet may be trapped i n a resona nce 


whic 


i can cause a significant increase in its orbital eccentricity flKlev fe Haghighipour 


20141 ). In the case of Kepler-34b the 10:1 mean motion resonace with the st ellar binary 


may; 


2015 ). 


have affected the evolution of the planetary eccentricity to some extent flChavez et al . 


7. DISCUSSION 

We would like to point out here that we assumed our six Kepler systems to be 
coplanar. In reality, they all have some mutual orbital inclination, with Kepler-413 having 
the largest one among those systems, i.e. / = 4.073°. In order to check whether our 
anayltic description of the system’s evolution remains valid for such inclinations, we ran 
two and three dimensional simulations for Kepler-413. We found that difference between 
the two models was insignificant for the amplitude of the eccentricity vector. Therefore, 
our predictions for the maximum and average eccentricities remain valid. However, the 
longitude of pericentre seemed to be affected considerably. Hence, we do not recommend to 
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m 0 


Fig. 1.— The Jacobi decomposition of the three body problem. 


Table 1: Masses and orbital elements for Kepler-16, Kepler-34, Kepler-35, Kepler-38, Kepler- 


64 and Kepler-413. 


System 

m o (M 0 ) 

mi(M 0 ) 

m 2 (Mj) 

ai (au) 

a 2 (au) 

ei 

Kepler-16 

o.6897ig;88ii 

0 909^+°-00066 
U.zuzdd_Q .00065 

fl ooq+0.016 

U.OOO_ooi6 

0 99/101+0.00035 

U.ZZ z ±Ol_o 00034 

0.70481°;°°" 

0 1 ^Q44+ 0 - 00061 

U.iOy^_ 0 .00062 

Kepler-34 

1 - 047918:8830 

1.0208™2 

0 220 +0 ' 011 
u.zzu_ 0 . 0 io 

O 99889+0-00019 

u.zzooz_ 0 00018 

1 . 089618 : 88 !! 

0 P )2D87 +0 ‘ 00052 

U.OZUOf _o.00055 

Kepler-35 

n 887£i+0-005i 
U.oo * ^—0.0053 

0.80941°;“" 

n 1 97+0-020 

U - 1Z ‘ -0.021 

n 1 1 7+0.00028 

-0.00029 

0 fiOS4 , j +0 ' 00100 

U.OU9 z ±O_ 0 .00102 

0 1491+ 0 - 0014 

U.l z ±Zi_o 0014 

Kepler-38 

0.949 

0.249 

<0.384 (95% conf.) 

0.1469 

0.4644 

0.1032 

Kepler-64 

1-38418:8?! 

0.3861°;"* 

<0.532 (99.7% conf.) 

n 1744+0.0031 

U.l ( ^_0.0031 

0.6341°;°" 

O 91 1 7+0-0051 
U.zil f —0 0051 

Kepler-413 

0.820 

0.5423 

0.21 

0.10148 

0.3553 

0.0365 
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Fig. 2.— Top row: logarithmic error (in color) for the averaged square planetary eccentricity 
(left) and the maximum planetary eccentricity (right). Bottom row: the logarithms of the 
numerical values (in colour) of the averaged squared planetary eccentricity (left) and the 
maximum planetary eccentricity (right). The mass parameters of the system are mi/(mo + 
rrii) = 0.5 and m/(m n + mi) = 10" 2 , X is the initial period ratio and X crit is the critical 


period ratio based on 
orbital period. 
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Fig. 3.— Top row: logarithmic error (in color) for the averaged square planetary eccentricity 
(left) and the maximum planetary eccentricity (right). Bottom row: the logarithms of the 
numerical values (in colour) of the averaged squared planetary eccentricity (left) and the 
maximum planetary eccentricity (right). The mass parameters of the system are mi /(mo + 
mi) = 0.3 and + m-i) = 10~ 3 , X is the initial period ratio and X crit is the critical 


period ratio based on 
secular period. 
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Fig. 4.— A zoom in of the plots of Figured for 1 < X/X crit < 2. 
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Fig. 5.— Top row: logarithmic error (in color) for the averaged square planetary eccentricity 
(left) and the maximum planetary eccentricity (right). Bottom row: the logarithms of the 
numerical values (in colour) of the averaged squared planetary eccentricity (left) and the 
maximum planetary eccentricity (right). The mass parameters of the system are mi /(mo + 
mi) = 0.1 and + m-i) = 10" 6 , X is the initial period ratio and X crit is the critical 


period ratio based on 
secular period. 
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time [secular periods] time [secular periods] 


Fig. 6.— Semi-major axis logarithmic relative change (left graph) and logarithmic absolute 
change in the eccentricity (right graph). The red colour represents the stellar orbit, while 
the green colour represents the planetary orbit. The mass parameters of the system are 
m 1 /(m 0 + mi) — 0.1 and + mi) = 10~ 2 , the initial stellar eccentricity is 0.5, the 

initial stellar eccentric anomaly is 0°, the initial planetary mean anomaly is 90°, the initial 
period ratio is 1.5 X cr i t and the integration time is one analytical secular period. Variations 
in the planetary eccentricity remain dominant even for systems close to the critical period 
ratio for dynamical instability. 
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Fig. 7.— Eccentricity against time for for a system with mo = 5M 0 , mi = 4M 0 , m 2 = 1 Mj, 
a>i = 0.2 au, a 2 = 6.84014558 an, ei = 0.5, W\ = 0°, E w = 0° and Z 2 0 = 90°. The red curve 
comes from the integration of the non-relativistic full equations of motion, the green curve 
comes from the integration of the relativistic full equations of motion and the blue curve is 
our analytic estimate. 
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Fig. 8.— Eccentricity against time for Kepler-16b, Kepler-34b and Kepler-35b. The red 
curve comes from the numerical integration of the full equations of motion, the green curve 
is our analytical estimates, the blue curve is the analytical secular solution, while the black 
line denotes the current value of the planetary eccentricity. The integration time is one 
planetary period for the left column and one analytical secular period for the right column. 
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Fig. 9.— Same as Figure 0 but for Kepler-38b, Kepler-64b and Kepler-413b. 
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Fig. 10.— Difference in the eccentricity (left) and in the longitude of pericentre (right) 
between a 2D and a 3D model of Kepler-413. See Figure [9] for a comparison with the current 
and predicted eccentricity estimates. The secular period for a coplanar planetary orbit is 
about 12 years. 
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use the fully analytic description presented in section [5] to describe a system’s long term 
orbit evolution for mutual inclinations around / ^ 5° or higher. Figure 9 presents the 
relevant information. 


8. SUMMARY 

In this work, we derived analytical estimates for the motion of coplanar circumbinary 
planets. The derivation was done for initially circular orbits and it assumed that the 
planetary eccentricity remains low (e .2 < 0.2). Short term motion was described by solving 
differential equations derived from the Runge-Lenz vector. Similarly, canonical perturbation 
theory was used to model the secular motion. Wokring with the eccentric instead of the 
true anomaly to describe the short term motion of the perturber leads to a solution that 
does not depend on a series expansion in terms of the eccentricity of the binary. This is 
a great advantage as we do not need to worry about the convergence of the solution and 
the magnitude of the eccentricity of the perturber. Hence, our derivations are valid for all 
eccentricities of the stellar binary as long as the planet is far enough from the binary so 
that a hierarchical triple system can be considered. 

In section [5] we have presented a complete analytical framework to describe the motion 
of a stellar binary and its circumbinary planet. Analytical expressions were derived for the 
averaged square eccentricity as well as for the maximum value of the planet’s eccentricity. 
The estimates were tested numerically for different stellar and planetary masses over a wide 
range of eccenricity and period ratios. The results showed very good agreement between the 
numerical simulations and the analytical estimates. Furthermore, in case it is required, the 
analytical formulae can easily be adjusted to incorporate some general relativistic effects. 
Finally, the analytical estimates were applied to six Kepler circumbinary systems with 
good results, as they provided valuable information about the evolution and possibly the 
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formation of those systems. 

The analytical estimates obtained in this work can be used, for instance, in analytic 
planet formation theories, in modelling transits in circumbinary planetary systems, or 
determing habitable zones around stellar binaries. 

S.E. is grateful for the support by the European Union Seventh Framework Program 
(FP7/2007-2013) under grant agreement no. 282703, as well as the travel funding by the 
Austrian FWF project S11608-N16 (sub-project of the NFN S116). 


APPENDIX 


The expressions for Pi(t) and P 2 {t) in equations (j3Tj) are presented. They contain the 
short periodic terms of the outer eccentric vector. For t — t 0 we have E\ = E w and l 2 = Z 20 - 
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NOTATION 

semimajor axis of the inner {i = 1) and outer (i = 2) orbit 
vacuum light speed 
integration constants 
eccentricity vectors 

x (j = 1) and y (j = 2) components of the inner and outer eccentricity vector 
eccentric anomalies 
true anomalies 
gravitational constant 

specific angular momentum of the outer orbit 
Hamiltonians 

mutual inclination of the binary and planetary orbits 
parameters of the long term solution 
Delaunnay actions 

Delaunnay angles (mean anomaly, argument of pericenter) 

masses of the two stars 

mass of the planet 

total mass of the system 

mass parameters 

mean motions 

Legendre polynomials 

x and y component terms of short period e 2 solution 

barycentric position vectors of the stars (i = 0,1) and the planet (i = 2) 

Jacobi vector of the stellar orbit 

Jacobi vector of the planetary orbit 

angle between r and R 

time 

stellar and planetary longitude of pericentre 

pericentre rotation matrices 

period ratio between outer and inner orbit 
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